home *** CD-ROM | disk | FTP | other *** search
-
-
-
- CCCCGGGGGGGGSSSSVVVVDDDD((((3333SSSS)))) CCCCGGGGGGGGSSSSVVVVDDDD((((3333SSSS))))
-
-
-
- NNNNAAAAMMMMEEEE
- CGGSVD - compute the generalized singular value decomposition (GSVD) of
- an M-by-N complex matrix A and P-by-N complex matrix B
-
- SSSSYYYYNNNNOOOOPPPPSSSSIIIISSSS
- SUBROUTINE CGGSVD( JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B, LDB,
- ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, RWORK,
- IWORK, INFO )
-
- CHARACTER JOBQ, JOBU, JOBV
-
- INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P
-
- INTEGER IWORK( * )
-
- REAL ALPHA( * ), BETA( * ), RWORK( * )
-
- COMPLEX A( LDA, * ), B( LDB, * ), Q( LDQ, * ), U( LDU, * ), V(
- LDV, * ), WORK( * )
-
- IIIIMMMMPPPPLLLLEEEEMMMMEEEENNNNTTTTAAAATTTTIIIIOOOONNNN
- These routines are part of the SCSL Scientific Library and can be loaded
- using either the -lscs or the -lscs_mp option. The -lscs_mp option
- directs the linker to use the multi-processor version of the library.
-
- When linking to SCSL with -lscs or -lscs_mp, the default integer size is
- 4 bytes (32 bits). Another version of SCSL is available in which integers
- are 8 bytes (64 bits). This version allows the user access to larger
- memory sizes and helps when porting legacy Cray codes. It can be loaded
- by using the -lscs_i8 option or the -lscs_i8_mp option. A program may use
- only one of the two versions; 4-byte integer and 8-byte integer library
- calls cannot be mixed.
-
- PPPPUUUURRRRPPPPOOOOSSSSEEEE
- CGGSVD computes the generalized singular value decomposition (GSVD) of an
- M-by-N complex matrix A and P-by-N complex matrix B:
- U'*A*Q = D1*( 0 R ), V'*B*Q = D2*( 0 R )
-
- where U, V and Q are unitary matrices, and Z' means the conjugate
- transpose of Z. Let K+L = the effective numerical rank of the matrix
- (A',B')', then R is a (K+L)-by-(K+L) nonsingular upper triangular matrix,
- D1 and D2 are M-by-(K+L) and P-by-(K+L) "diagonal" matrices and of the
- following structures, respectively:
-
- If M-K-L >= 0,
-
- K L
- D1 = K ( I 0 )
- L ( 0 C )
- M-K-L ( 0 0 )
-
- K L
-
-
-
- PPPPaaaaggggeeee 1111
-
-
-
-
-
-
- CCCCGGGGGGGGSSSSVVVVDDDD((((3333SSSS)))) CCCCGGGGGGGGSSSSVVVVDDDD((((3333SSSS))))
-
-
-
- D2 = L ( 0 S )
- P-L ( 0 0 )
-
- N-K-L K L
- ( 0 R ) = K ( 0 R11 R12 )
- L ( 0 0 R22 )
- where
-
- C = diag( ALPHA(K+1), ... , ALPHA(K+L) ),
- S = diag( BETA(K+1), ... , BETA(K+L) ),
- C**2 + S**2 = I.
-
- R is stored in A(1:K+L,N-K-L+1:N) on exit.
-
- If M-K-L < 0,
-
- K M-K K+L-M
- D1 = K ( I 0 0 )
- M-K ( 0 C 0 )
-
- K M-K K+L-M
- D2 = M-K ( 0 S 0 )
- K+L-M ( 0 0 I )
- P-L ( 0 0 0 )
-
- N-K-L K M-K K+L-M
- ( 0 R ) = K ( 0 R11 R12 R13 )
- M-K ( 0 0 R22 R23 )
- K+L-M ( 0 0 0 R33 )
-
- where
-
- C = diag( ALPHA(K+1), ... , ALPHA(M) ),
- S = diag( BETA(K+1), ... , BETA(M) ),
- C**2 + S**2 = I.
-
- (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N), and R33 is stored
- ( 0 R22 R23 )
- in B(M-K+1:L,N+M-K-L+1:N) on exit.
-
- The routine computes C, S, R, and optionally the unitary
- transformation matrices U, V and Q.
-
- In particular, if B is an N-by-N nonsingular matrix, then the GSVD of A
- and B implicitly gives the SVD of A*inv(B):
- A*inv(B) = U*(D1*inv(D2))*V'.
- If ( A',B')' has orthnormal columns, then the GSVD of A and B is also
- equal to the CS decomposition of A and B. Furthermore, the GSVD can be
- used to derive the solution of the eigenvalue problem:
- A'*A x = lambda* B'*B x.
- In some literature, the GSVD of A and B is presented in the form
- U'*A*X = ( 0 D1 ), V'*B*X = ( 0 D2 )
-
-
-
- PPPPaaaaggggeeee 2222
-
-
-
-
-
-
- CCCCGGGGGGGGSSSSVVVVDDDD((((3333SSSS)))) CCCCGGGGGGGGSSSSVVVVDDDD((((3333SSSS))))
-
-
-
- where U and V are orthogonal and X is nonsingular, and D1 and D2 are
- ``diagonal''. The former GSVD form can be converted to the latter form
- by taking the nonsingular matrix X as
-
- X = Q*( I 0 )
- ( 0 inv(R) )
-
-
- AAAARRRRGGGGUUUUMMMMEEEENNNNTTTTSSSS
- JOBU (input) CHARACTER*1
- = 'U': Unitary matrix U is computed;
- = 'N': U is not computed.
-
- JOBV (input) CHARACTER*1
- = 'V': Unitary matrix V is computed;
- = 'N': V is not computed.
-
- JOBQ (input) CHARACTER*1
- = 'Q': Unitary matrix Q is computed;
- = 'N': Q is not computed.
-
- M (input) INTEGER
- The number of rows of the matrix A. M >= 0.
-
- N (input) INTEGER
- The number of columns of the matrices A and B. N >= 0.
-
- P (input) INTEGER
- The number of rows of the matrix B. P >= 0.
-
- K (output) INTEGER
- L (output) INTEGER On exit, K and L specify the dimension
- of the subblocks described in Purpose. K + L = effective
- numerical rank of (A',B')'.
-
- A (input/output) COMPLEX array, dimension (LDA,N)
- On entry, the M-by-N matrix A. On exit, A contains the
- triangular matrix R, or part of R. See Purpose for details.
-
- LDA (input) INTEGER
- The leading dimension of the array A. LDA >= max(1,M).
-
- B (input/output) COMPLEX array, dimension (LDB,N)
- On entry, the P-by-N matrix B. On exit, B contains part of the
- triangular matrix R if M-K-L < 0. See Purpose for details.
-
- LDB (input) INTEGER
- The leading dimension of the array B. LDB >= max(1,P).
-
- ALPHA (output) REAL array, dimension (N)
- BETA (output) REAL array, dimension (N) On exit, ALPHA and
- BETA contain the generalized singular value pairs of A and B;
-
-
-
- PPPPaaaaggggeeee 3333
-
-
-
-
-
-
- CCCCGGGGGGGGSSSSVVVVDDDD((((3333SSSS)))) CCCCGGGGGGGGSSSSVVVVDDDD((((3333SSSS))))
-
-
-
- ALPHA(1:K) = 1,
- BETA(1:K) = 0, and if M-K-L >= 0, ALPHA(K+1:K+L) = C,
- BETA(K+1:K+L) = S, or if M-K-L < 0, ALPHA(K+1:M)= C,
- ALPHA(M+1:K+L)= 0
- BETA(K+1:M) = S, BETA(M+1:K+L) = 1 and ALPHA(K+L+1:N) = 0
- BETA(K+L+1:N) = 0
-
- U (output) COMPLEX array, dimension (LDU,M)
- If JOBU = 'U', U contains the M-by-M unitary matrix U. If JOBU =
- 'N', U is not referenced.
-
- LDU (input) INTEGER
- The leading dimension of the array U. LDU >= max(1,M) if JOBU =
- 'U'; LDU >= 1 otherwise.
-
- V (output) COMPLEX array, dimension (LDV,P)
- If JOBV = 'V', V contains the P-by-P unitary matrix V. If JOBV =
- 'N', V is not referenced.
-
- LDV (input) INTEGER
- The leading dimension of the array V. LDV >= max(1,P) if JOBV =
- 'V'; LDV >= 1 otherwise.
-
- Q (output) COMPLEX array, dimension (LDQ,N)
- If JOBQ = 'Q', Q contains the N-by-N unitary matrix Q. If JOBQ =
- 'N', Q is not referenced.
-
- LDQ (input) INTEGER
- The leading dimension of the array Q. LDQ >= max(1,N) if JOBQ =
- 'Q'; LDQ >= 1 otherwise.
-
- WORK (workspace) COMPLEX array, dimension (max(3*N,M,P)+N)
-
- RWORK (workspace) REAL array, dimension (2*N)
-
- IWORK (workspace/output) INTEGER array, dimension (N)
- On exit, IWORK stores the sorting information. More precisely,
- the following loop will sort ALPHA for I = K+1, min(M,K+L) swap
- ALPHA(I) and ALPHA(IWORK(I)) endfor such that ALPHA(1) >=
- ALPHA(2) >= ... >= ALPHA(N).
-
- INFO (output) INTEGER
- = 0: successful exit.
- < 0: if INFO = -i, the i-th argument had an illegal value.
- > 0: if INFO = 1, the Jacobi-type procedure failed to converge.
- For further details, see subroutine CTGSJA.
-
- PPPPAAAARRRRAAAAMMMMEEEETTTTEEEERRRRSSSS
- TOLA REAL
- TOLB REAL TOLA and TOLB are the thresholds to determine the
- effective rank of (A',B')'. Generally, they are set to TOLA =
- MAX(M,N)*norm(A)*MACHEPS, TOLB = MAX(P,N)*norm(B)*MACHEPS. The
-
-
-
- PPPPaaaaggggeeee 4444
-
-
-
-
-
-
- CCCCGGGGGGGGSSSSVVVVDDDD((((3333SSSS)))) CCCCGGGGGGGGSSSSVVVVDDDD((((3333SSSS))))
-
-
-
- size of TOLA and TOLB may affect the size of backward errors of
- the decomposition.
-
- Further Details ===============
-
- 2-96 Based on modifications by Ming Gu and Huan Ren, Computer
- Science Division, University of California at Berkeley, USA
-
- SSSSEEEEEEEE AAAALLLLSSSSOOOO
- INTRO_LAPACK(3S), INTRO_SCSL(3S)
-
- This man page is available only online.
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- PPPPaaaaggggeeee 5555
-
-
-
-
